Differential Equations

Sesame, Door Opens.

As well-known history, Sir Newton used the tools of Calculus to study mechanics.

Colan

Only one Truth there

Sympy: Symbolic Differential Equation Solver

Newton's Cooling Law

Let \begin{eqnarray} T_0 & : & \text{the initial temperature of a body at } t_0\ T (t) & : & \text{the temperature of the body at time $t$,}\ T_m & : & \text{the temperature of certain a medium and $T_m < T_0$.} \end{eqnarray} Now put the body into the medium. The cooling law states that the rate which the body cools at any given instant is proportional to the difference between its at that instant and the the temperature $T_m$ of the surrounding medium, i.e. $$ \frac{T(t+\Delta t)-T(t)}{\Delta t}\sim T(t)-T_m $$

Note

\begin{eqnarray} & T' (t) & = k (T (t) - T_m)\ & T (t_0) & = T_0\ \Longrightarrow & \int \frac{d T}{T - T_m} & = \int k d t\ \Longrightarrow & \ln (T - T_m) & = k t + C\ \Longrightarrow & T & = T_m + c e^{k t}\ \Longrightarrow & T & = T_m + (T_0 - T_m) e^{ k (t-t_0)} \end{eqnarray} by the initial conditions.

In [1]:
from sympy import *
from sympy.abc import t,x,k,N

init_printing()
In [2]:
var('C1 C2 C3')
y=Function("y")
def Dsol(func,t0,initial):
    """
    ODE solver
    """
    y=Function("y")
    yp=Derivative(y(t),t)
    sol=dsolve(yp-func,y(t))
    cc=solve(sol.rhs.subs(t,t0)-initial,C1)
    return sol.rhs.subs(C1,cc[0])
In [3]:
var("L t0 T0 Tm")
func=k*(y(t)-Tm)

Dsol(func,t0,T0)
Out[3]:
$$Tm + \left(T_{0} - Tm\right) e^{k t} e^{- k t_{0}}$$

Question

A dead body was found around Carven Garden, London. The Holmes and Dr. Watson arrived the scene at 9 : 00 A.M. and immediately recorded the temperature of the body as $30.1^\circ$C. One hour later they noted that the body's temperature had dropped to $29.2^\circ$C, and during the hour, room temperature had remained at $20^\circ$ˆ˜C. Assume that the body's temperature obeys Newton's law of cooling, determine when death occurred. (Normal body temperature is $37^\circ$C.)

Answer:

Which is the following is right? a) AM 2hr and 09min b) AM 3hr and 26min c) AM 4hr and 47min

Pursiut Curves

The path an object takes when chasing after another object in the most effective way.

In Francis Godwin’s story, The Man in the Moone, an astronaut harnesses a wedge of 25 swans and flies to the moon. The swans fly at a constant rate and always head toward the moon, in accordance with their annual migration. Thus, the trajectory from the earth to the moon is not a straight line, but is a pursuit curve. In the story, Godwin says that the flight to the moon takes twelve days, whereas the return trip, which follows a straight line, takes eight days days.

How good a guess was Godwin’s?

References:

  1. Francis Godwin (1562~1633) , The Man in the Moone
  2. A. Simonson, Pursuit Cruves for the Man in the Moon, The Mathematical Association of America, (?)
In [ ]:

Problems Exploring and Showing

0. Basics of Numerical: from Continuous to discrete
1. Basics of Physics;
2. Pictures describing the tractories of observerved objects;
3. Animation of whole the processes;
4. Simple Interaction by passing the Parameter to System.

Skills Requirement

1. Python >3.x
2. numpy: data handling
3. matplotlib-1.3.1: making pictures and animation,
4. IPython, JSAnimation, IPyWidget: make results well display and manipulation from HTML-based browsers.
In [ ]:

Differential Equations

Equations involve differential operation. Nowaday, it is widely studied in theoretic and applied science.

How to Solve DE's

Not all differential equations could be solved. Could not solve it, then simulate it out!

Difference Scheme

Consider the following initial value problem: \begin{array}{cccl} &x'(t)&=&x(t), \cr &x(0)&=&1 \cr \Rightarrow &x &=& \exp(t). \end{array}
The respective numerical apprroximation, $$x'(t) \sim \frac {x(t+\vartriangle t)-x(t)}{\vartriangle t}.$$

Now partition $[0,t]$ into $n$ equally-like nonoverlapping sub-intervals as follows:

\begin{array}{cccccc} x: & x_0 = 1 & x_1 & x_2 & \cdots & x_{n-1} & \quad x_n \cr & \bot \ldots & ... \perp ... & ... \perp ... & \ldots \ldots . & ... \perp ... & \ldots \bot\cr t: & 0\quad & \vartriangle t & 2\vartriangle t & \cdots & (n-1)\vartriangle t & \quad n\vartriangle t=t \end{array}

where $x_i=x(i\vartriangle t),$ and $\vartriangle t=\frac{t}{n}$.

This induces the following difference scheme: \begin{array}{cccl} &x'&=& x \cr &\frac{x_{i+1}-x_i}{\vartriangle t} &\sim& xi \cr \Rightarrow &x{i+1} &\sim& (1+\vartriangle t)xi &\sim& (1+\vartriangle t)^2x{i-1} \ &&\sim& (1+\vartriangle t)^{i+1} x_0 \end{array} for $i=1,2,\cdots$.

$a^\circ$) Estimate the value of Euler number,$e\sim2.71828\cdots$
$e=\exp(1)$ could be evaluated by \begin{array}{ccl} x_0 &=& e(0), \cr x_1 &=& e(\vartriangle t) \cr &\sim& (1+\vartriangle t)x_0, \cr x_2 &=& e(2\vartriangle t) \cr &\sim& (1+\vartriangle t)x_1 \cr &\sim& (1+\vartriangle t)^2 x_0 \cr &\vdots \cr x_n &\sim& (1+\vartriangle t)^n x_0 \end{array}
, which refers to the following partition for $\exp(x)$ from 0 to 1:

\begin{array}{cccccc} x: & x_0 = 1 & x_1 & x2 & \cdots & x{n-1} & x_n = e\cr &\bot \ldots & ... \perp ... & ... \perp ... & \ldots \ldots . & ... \perp ... & \ldots \bot\cr t: & 0\quad & \frac{1}{n} & \frac{2}{n} & \cdots & \frac{n-1}{n} & \quad1 \end{array}

where $\vartriangle t=\frac{1}{n}$. And the approximation for Euler number, $e$, could be estimeted by the following: $$ e=x_n\sim \left(1+\frac{1}{n}\right)x_{n-1}\sim \left(1+\frac{1}{n}\right)^2x_{n-2}\sim\dots \sim \left(1+\frac{1}{n}\right)^n e(0)=\left(1+\frac{1}{n}\right)^n$$

In [4]:
n=100
a=0;b=1;
dt=(b-a)/float(n)
e_value=(1+dt)**n
print 'the approximation of e is ', e_value

the approximation of e is  2.70481382942
In [5]:
from IPython.display import clear_output
import sys
import time

def sim(n):
    a=0;b=1;
    dt=(b-a)/float(n)
    e_value=(1+dt)**n
    print "the approximation, ( 1  +  1/",n,")^",n,",  of e is " ,e_value

for i in range(1,11):
    time.sleep(0.5)
    clear_output()
    n=10*i
    sim(n)
    sys.stdout.flush()

the approximation, ( 1  +  1/ 100 )^ 100 ,  of e is  2.70481382942
In [ ]:
In [6]:
%matplotlib inline
import numpy as np
from numpy import exp
import matplotlib.pylab as plt

"""
from pylab import *
rcParams['figure.figsize'] = 4,4
"""

n=11
a=0.;b=1.
dx=(b-a)/float(n-1)
x=np.linspace(a,b,n)

y=np.ones(n)
for i in np.arange(n):
      y[1:]=(1+dx)*y[:-1]
P2=plt.plot(exp(x),'b-',y,'r*')

$b^\circ$) Logistic Model (or called gassip model, learning model,...)

\begin{align} y'&=&a y (M-y) \cr y(0) &=& y_0 \end{align}


where $a,M, y_0$ are constants and $M$ is the maximal resources that environment could afford.

The scheme runs as follows: \begin{array}{ccl} &\frac{y_{n+1}-y_n}{\vartriangle t} &=& a y_n (M-yn) \cr \Rightarrow &y{n+1} &=& y_n +a y_n (M-y_n)\vartriangle t \end{array}

The solution curve is like $S$-shape:

In [8]:
a=1/5.
M=20
t0=0;t1=1
n=30
u0=5
t=np.linspace(t0,t0+t1,n+1)
dt=(t1-t0)/float(n)

u=np.zeros(np.size(t))
u=u0*np.ones(np.size(u))
u_h=M/2.*np.ones(np.size(u))

for i in np.arange(len(u)-1):
    """
    i.e.
    u[i+1] = u[i]+a*u[i]*(M-u[i])*dt
    """
    u[1:]=u[:-1]+a*u[:-1]*(M-u[:-1])*dt

plt.plot(t,u,t,u_h,'r--')
Out[8]:
[<matplotlib.lines.Line2D at 0x1061e2c90>,
 <matplotlib.lines.Line2D at 0x1061e2f10>]
In [11]:
Levels=[1,2,5,10,13,15,17,25,22,21]
P=[]

for level in Levels:
    u=np.zeros(np.size(t))
    u[0]=level
    for i in np.arange(len(u)-1):
        u[1:]=u[:-1]+a*u[:-1]*(M-u[:-1])*dt
    P+=plt.plot(t,u)
P+=plt.plot(u_h,'r--')
P+=plt.plot(2*u_h,'b--')
plt.axis([0,1,0,26])
#show(P)
Out[11]:
[0, 1, 0, 26]

Question

Newton's law of cooling states that the rate change of temperate of heated object is propotional to the difference between its own temperature and the ambient temperature. $$ \mathbf{y(t)=k(y(t)-y_m)}$$ where $y(0)=y_0$ is initial temperature of temperature, $y'(t)$ is the temperature at time $t$, and $y_m$ is ambient temperature. Suppose that a corpse was discovered in a room and its temperature was 32°C. The temperature of the room is kept constant at 22°C. Three hours later the temperature of the corpse dropped to 27°C.

  1. Find the time of death.
  2. make a picture about the temperature curve of corpse change.

SciPy

Known as Python Repository library of Scientific Computations, it includes the schemes which are capable of doing numeric ODE's like above examples.

Example Why is there no people died of being hit by falling rain?

  • According Newton's Law, the rain drop will fall down as a free body and be formulated by the following: $$\ddot{\mathbf{x}} = \mathbf{g}-r|\dot{\mathbf{x}}|^2 $$

    • initial conditions, $(\mathbf{x},\dot{\mathbf{x}})=(\mathbf{x},\mathbf{v})=(0,0)$;
    • $\mathbf{g}=9.81 m/s^2$ is the gravity acceleration;
    • the air resistant force $r|\dot{\mathbf{x}}|^2$

Solved by SciPy

  • $$\left[\matrix{\dot x\cr \ddot x}\right]=\dot{\left[\matrix{x\cr v}\right]}=\frac{d}{ d t}\left[\matrix{x\cr v}\right]$$
  • $$\frac{d }{ d t}\mathbf{X}=\frac{d }{ d t}{\left[\matrix{x\cr v}\right]}={\left[\matrix{v\cr g-r v^2}\right]}=f(\mathbf{X})$$
  • Solved by scipy.integrate.odeint(f,x_init,t)
In [ ]:
In [1]:
import numpy as np
import scipy.integrate as sp
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
# The initial position and velocity is (0,0) (or * above sea level)
x0 = np.zeros(2)
In [3]:
g=9.81
r=1.
In [6]:
def f(x,t, r):
    # X has twocomponents: X=[x, x'].
    x, xdot = x[:1], x[1:]
    # We compute the second derivative x'' of x.
    xdotdot = g-r * xdot *xdot
    # We return X'=[x', x''].
    return np.r_[xdot, xdotdot]
In [7]:
t = np.linspace(0., 10., 100)
x = sp.odeint(f, x0, t, args=(r,))
In [8]:
from numpy import sqrt
plt.plot(t,sqrt(g)*np.ones(np.size(t)))
plt.plot(t[::2],x[::2,1], 'o-', mew=1, ms=8, mec='w',mfc="r")
Out[8]:
[<matplotlib.lines.Line2D at 0x1064c7c90>]

Conclusion and Discussion

  • the velocity, $v=\dot x$, seems to approach about $\mathbf{\sqrt g}$ m/sec here ( at which $\frac{d }{d t}\dot v=0$).
  • Actually, we also have
    • $v=\frac{ d x}{d t}$
    • $$ \mathbf{\dot v}=\mathbf{g} - r \mathbf{v^2}$$
      • $\dot v >0$ if $ g-r v^2>0$ (i.e. $v<\sqrt{g/r}$ $\Longrightarrow v$ is increasing;
      • $\dot v <0$ if $ g-r v^2<0$ (i.e. $v>\sqrt{g/r}$ $\Longrightarrow v$ is decreasing;
      • $v$ is increasing from initial 0 under the condition $g-r v^2: >0$ to $0$
    • The exact solution of $v$: $$ v= \frac{2\sqrt{g}}{1+\exp(-2\sqrt g t)}-\sqrt g$$ and $$ \lim_{t\to\infty}v(t)=\sqrt g$$
    • Quantitative method: easy to use to make analysis
      • stable point: $\mathbf{g} - r \mathbf{v^2}=0\to v=\sqrt {g/r}$;
      • $\dot v=g-rv^2 >0$ if $v<\sqrt{g/r}$ $\Longrightarrow v$ is increasing;
      • $\ddot v =-2rv<0 => v$ is concave downwards.
In [28]:
from sympy import Symbol , dsolve , Function , Derivative ,symbols, Symbol,solve, simplify,init_printing,integrate,exp
#from sympy.abc import a,r,t
init_printing()
In [37]:
t = Symbol("t",positive=True)
r,x = symbols("r x",positive=True)
v=Function('v')
g=9.8
#r=1
f_=Derivative(v(t),t)-g+r*v(t);print f_

r*v(t) + Derivative(v(t), t) - 9.8
In [40]:
sol=dsolve(f_,v(t));print sol.rhs

(1.0*exp(r*(C1 - t)) + 9.8)/r
In [ ]:

Dog Tracking

Suppose that a dog, moving at constant speed while always pointing towards its master.

Consider:

  • $\vec{v_D}$: constant velocity of dog, 10 m/s

  • $\vec{v_M}$: constant velocity of master, 5 m/s, and towards east.

    </ul> </font> Initially, the master stands at the north of dog and 10 meters away it.

    The motion could be represented simply by the vector form: $$\vec{v_D}=|\vec{v_D}| \cdot \frac{\vec{M}-\vec{D}}{||\vec M-\vec D||}$$
    Then the coordinates of position is as follows:

    $$ \vec M(\vec D)(t_{n+1})=\vec M(\vec D)(t_n)+vel_{M(D)}(t_n)\times d t$$

Outline

  • Dog positions: Dpos = pos1, pos2, ...
  • Master position: Mpos = pos1, pos2 ,...
In [1]:
In [2]:
%matplotlib inline

import numpy as np
from numpy import pi,sin,cos,size
from numpy.linalg import norm
import matplotlib.pyplot as plt
from matplotlib import animation
from JSAnimation import IPython_display
from ipywidgets import interact, IntSlider#,Radio
from IPython.display import display, clear_output
In [2]:
Dpos=np.array([[0,0]])
Mpos=np.array([[0,10]])
In [3]:
dt=0.1
Mvel=np.array([5,0])
In [4]:
#while i<40:
#while norm(Dpos[-1]-Mpos[-1])>0.01:
while Dpos[-1][1]<Mpos[0][1]:
  Dvel=10*(Mpos[-1]-Dpos[-1])/norm(Mpos[-1]-Dpos[-1])
  
  Mpos1 = Mpos[-1]+Mvel*dt
  Dpos2 = Dpos[-1]+Dvel*dt


  Mpos=np.vstack([Mpos,Mpos1])
  Dpos=np.vstack([Dpos,Dpos2]) 
  #i=i+1
In [5]:
figM = plt.figure(figsize=(6,6))
ax = figM.add_subplot(111)
ax.set_title("Dog Chasing",fontsize=14)

plt.xlim([0,12])
plt.ylim([0,12])
ax.scatter(Mpos[:,0],Mpos[:,1],  c='r' )
ax.scatter(Dpos[:,0],Dpos[:,1],c='g')
Out[5]:
<matplotlib.collections.PathCollection at 0x107f6b438>
In [8]:
def DogTrack(Vel):
    #while i<40:
    #while norm(Dpos[-1]-Mpos[-1])>0.01:
    Dpos=np.array([[0,0]])
    Mpos=np.array([[0,10]])
    i=1
    #while Dpos[-1][1]<Mpos[0][1] or i<40:
    while i<40:    
      Dvel=Vel*(Mpos[-1]-Dpos[-1])/norm(Mpos[-1]-Dpos[-1])
  
      Mpos1 = Mpos[-1]+Mvel*dt
      Dpos2 = Dpos[-1]+Dvel*dt

      Mpos=np.vstack([Mpos,Mpos1])
      Dpos=np.vstack([Dpos,Dpos2]) 
      i=i+1
          
    figM = plt.figure(figsize=(9,7))
    #ax = figM.add_subplot(211)
    plt.title("Dog Chasing",size=14)

    plt.xlim([0,18])
    plt.ylim([0,14])
    plt.scatter(Mpos[:,0],Mpos[:,1], c='r' )
    plt.scatter(Dpos[:,0],Dpos[:,1],c='g')
    plt.text(1, 11.5,
            "Master Velocity={0}\nDog Velocity={1:.2f}".format(5,Vel),
            size=14, color='gray')
    #return figM
In [9]:
interact(DogTrack,Vel=IntSlider(min=2, max=10, step=1))
In [ ]:

Cats Chasing

Four cats are in the four corners a squareroom of unitary base. At time t=0, they are starting to chase the tails of others. Same as the dog tracking problem:

The case of the cat at the left lower corner runs twice faster than others is simulated as follows:

$$vel_i=\frac{\vec{P}_i-\vec{P}_{i-1}}{||\vec{P}_i-\vec{P}_{i-1}||}$$$$ \vec P_{i}(t_{n+1})=\vec P_{i}(t_n)+vel_{i}(t_n)\times d t$$
In [3]:
Cpos1=np.array([[0,0]])
Cpos2=np.array([[1,0]])
Cpos3=np.array([[1,1]])
Cpos4=np.array([[0,1]])
In [4]:
dt=0.1
i=1
In [5]:
#while i<40:
while norm(Cpos1[-1]-Cpos2[-1])>0.01:
  v1=2*(Cpos2[-1]-Cpos1[-1])
  v2=2*(Cpos3[-1]-Cpos2[-1])
  v3=2*(Cpos4[-1]-Cpos3[-1])
  v4=3*(Cpos1[-1]-Cpos4[-1])

  Cpos11 = Cpos1[-1]+v1*dt
  Cpos21 = Cpos2[-1]+v2*dt
  Cpos31 = Cpos3[-1]+v3*dt
  Cpos41 = Cpos4[-1]+v4*dt 

  Cpos1=np.vstack([Cpos1,Cpos11])
  Cpos2=np.vstack([Cpos2,Cpos21]) 
  Cpos3=np.vstack([Cpos3,Cpos31])
  Cpos4=np.vstack([Cpos4,Cpos41])
In [6]:
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
ax.set_title("Four Cats Pursuit Curve",fontsize=14)

plt.xlim([0,1])
plt.ylim([0,1])
ax.scatter(Cpos1[:,1],Cpos1[:,0], c='y' )
ax.scatter(Cpos2[:,1],Cpos2[:,0],c='g')
ax.scatter(Cpos3[:,1],Cpos3[:,0],c='b')
ax.scatter(Cpos4[:,1],Cpos4[:,0],c='r')
Out[6]:
<matplotlib.collections.PathCollection at 0x110036668>
In [7]:
anim = plt.figure(figsize=(6,6))
ax = anim.add_subplot(111)
ax.set_title("Four Cats Pursuit Curve",fontsize=14)

plt.xlim([0,1])
plt.ylim([0,1])

def init():
    return ax.scatter(Cpos1[0,1],Cpos1[0,0], c='y' ), \
           ax.scatter(Cpos2[0,1],Cpos2[0,0],c='g'), \
           ax.scatter(Cpos3[0,1],Cpos3[0,0],c='b'), \
           ax.scatter(Cpos4[0,1],Cpos4[0,0],c='r')
def animate(i):
    return ax.scatter(Cpos1[0:i,1],Cpos1[0:i,0], c='y' ), \
       ax.scatter(Cpos2[0:i,1],Cpos2[0:i,0],c='g'), \
       ax.scatter(Cpos3[0:i,1],Cpos3[0:i,0],c='b'), \
       ax.scatter(Cpos4[0:i,1],Cpos4[0:i,0],c='r')
 
animation.FuncAnimation(anim, animate, init_func=init, frames=size(Cpos1[:,1]), interval=size(Cpos1[:,1]))
Out[7]:


Once Loop Reflect

One run faster

Suppose that one is runs faster than others, Then ... This is nothing to be worry because of the trusted simulation scheme doing well.

In [15]:
Cpos1=np.array([[0,0]])
Cpos2=np.array([[1,0]])
Cpos3=np.array([[1,1]])
Cpos4=np.array([[0,1]])

#while i<40:
while norm(Cpos1[-1]-Cpos2[-1])>0.01:
  v1=(Cpos2[-1]-Cpos1[-1])
  v2=(Cpos3[-1]-Cpos2[-1])
  v3=(Cpos4[-1]-Cpos3[-1])
  v4=(Cpos1[-1]-Cpos4[-1])

  Cpos11 = Cpos1[-1]+3.*v1*dt
  Cpos21 = Cpos2[-1]+v2*dt
  Cpos31 = Cpos3[-1]+v3*dt
  Cpos41 = Cpos4[-1]+v4*dt 

  Cpos1=np.vstack([Cpos1,Cpos11])
  Cpos2=np.vstack([Cpos2,Cpos21]) 
  Cpos3=np.vstack([Cpos3,Cpos31])
  Cpos4=np.vstack([Cpos4,Cpos41])
In [16]:
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
ax.set_title("Four Cats Pursuit Curve",fontsize=14)

plt.xlim([0,1])
plt.ylim([0,1])
ax.scatter(Cpos1[:,1],Cpos1[:,0], c='y' )
ax.scatter(Cpos2[:,1],Cpos2[:,0],c='g')
ax.scatter(Cpos3[:,1],Cpos3[:,0],c='b')
ax.scatter(Cpos4[:,1],Cpos4[:,0],c='r')
Out[16]:
<matplotlib.collections.PathCollection at 0x1090e44e0>

Exercise:

  1. Suppose that there are three cats stand at the corners of the equilateral triangle. Guess the trajectories about the cats chansing others' tails.

  2. If they run wih the different speed with respect to others, what would happen?

In [ ]:

Travel to the Moon

In 1968, Apollo 8, the first manned craft to circle the moon, reached the moon in 3.5 days, and returned in 2.5 days.

In [ ]:

Data

How fast does the Moon travel around Earth?

The Moon orbits Earth at a speed of 2,288 miles per hour (3,683 kilometers per hour). During this time it travels a distance of 1,423,000 miles (2,290,000 kilometers).

$$ 1 \text{ Lunar distance }=1 LU\sim 384,000 \text{ km }\\ p \text{ period } \sim 29.5 \text{ days}\\ R_{\text{earth}}~0.067 LU, r_{\text{moon}}~0.00453 LU \\ M(t): \text{ position of Moon }= (M_1(t),M_2(t))=\left(\cos\left(\frac{2\pi t}{p}\right), \sin\left(\frac{2\pi t}{p}\right)\right)\\ (x,y): \text{ position of spacecraft} $$

Differential Equations:

$$ \frac{y'}{x'}=\frac{M_2-y}{M_1-x}\\ (x')^2+(y')^2=c^2 \text{ where c: velocity of spacecraft}$$

where $(x,y)$ is coordinate of spacecraft and $(M_1,M_2)$ is the position of Lunar.

In [17]:
Mpos=np.array([[1,0]])
Epos=np.array([[0,0]])
p=29.5
c=0.15
dt=0.1
i=1
In [18]:
while i<400:
#while norm(Dpos[-1]-Mpos[-1])>0.01:
#while Dpos[-1][1]<Mpos[0][1]:
  Mpos1=np.array([cos(2*pi*i*dt/p),sin(2*pi*i*dt/p)])  
  Evel=c*(Mpos[-1]-Epos[-1])/norm(Mpos[-1]-Epos[-1])
  #print norm(Mpos[-1]-Epos[-1])
  Epos2 = Epos[-1]+Evel*dt
  #print Evel,Epos2
  Mpos=np.vstack([Mpos,Mpos1])
  Epos=np.vstack([Epos,Epos2]) 
  i=i+1
In [19]:
figM = plt.figure(figsize=(6,6))
ax = figM.add_subplot(111)
ax.set_title("Moon Landing",fontsize=14)

plt.xlim([-1.1,1.1])
plt.ylim([-1.1,1.1])
ax.scatter(Mpos[:,0],Mpos[:,1], c='r' )
ax.scatter(Mpos[-1,0],Mpos[-1,1], c='b',s=200)
ax.scatter(Epos[-1,0],Epos[-1,1], c='b',s=80)
ax.plot([Mpos[-1,0],Epos[-1,0]],[Mpos[-1,1],Epos[-1,1]],'b-->')
#ax.quiver(Mpos[-1,0]-0.2, Mpos[-1,1]-0.2, (Mpos[-1,0]-Epos[-1,0]), (Mpos[-1,1]-Epos[-1,1]), pivot='middle',\
#          headwidth=1, headlength=1,color="red")

ax.scatter(Epos[:,0],Epos[:,1],c='g')
Out[19]:
<matplotlib.collections.PathCollection at 0x1093947b8>
In [20]:
Epos[-1,0], Epos[-1,1], Mpos[-1,0]-Epos[-1,0], Mpos[-1,1]-Epos[-1,1]
Out[20]:
(0.10286913269798306,
 0.69496494915895923,
 -0.70350224211884416,
 0.10455982710968648)
In [21]:
anim = plt.figure(figsize=(6,6))
ax = anim.add_subplot(111)
ax.set_title("Moon Landing",fontsize=14)
ax.text(-0.45, -0.5,
            "Velocity of Spacecraft={0:.2f}".format(c),
            fontsize=14, color='gray')

plt.xlim([-1.1,1.1])
plt.ylim([-1.1,1.1])

def init():
    return ax.scatter(Mpos[0,0],Mpos[0,1], c='r' ), \
       ax.scatter(Epos[0,0],Epos[0,1],c='g')
def animate(i):
    return ax.scatter(Mpos[0:4*i:4,0],Mpos[0:4*i:4,1], c='r' ), \
       ax.scatter(Epos[0:4*i:4,0],Epos[0:4*i:4,1],c='g'), \
       ax.plot([Mpos[4*i,0],Epos[4*i,0]],[Mpos[4*i,1],Epos[4*i,1]],'b--')

    
animation.FuncAnimation(anim, animate, init_func=init, frames=size(Mpos[::4,1]), interval=size(Mpos[::4,1]))
Out[21]:


Once Loop Reflect
In [27]:
def MoonTraj(vec):
    c=vec
    Mpos=np.array([[1,0]])
    Epos=np.array([[0,0]])

    i=1
    #while Dpos[-1][1]<Mpos[0][1] or i<40:
    while i<800:
          Mpos1=np.array([cos(2*pi*i*dt/p),sin(2*pi*i*dt/p)])  
          Evel=c*(Mpos[-1]-Epos[-1])/norm(Mpos[-1]-Epos[-1])

          Epos2 = Epos[-1]+Evel*dt

          Mpos=np.vstack([Mpos,Mpos1])
          Epos=np.vstack([Epos,Epos2]) 
          i=i+1

    figM = plt.figure(figsize=(6,6))
    #ax = figM.add_subplot(111)
    plt.title("Moon Landing",size=14)

    plt.xlim([-1.1,1.1])
    plt.ylim([-1.1,1.1])
    plt.scatter(Mpos[:,0],Mpos[:,1], c='r' )
    plt.scatter(Epos[::4,0],Epos[::4,1],c='g')            
    plt.text(-0.5, -0.5,
            "Velocity of Spacecraft={0:.2f}".format(c),
            size=14, color='gray')       
            
    #return figM
In [28]:
from ipywidgets import FloatSlider
interact(MoonTraj,vec=FloatSlider(min=0.05, max=0.30, step=0.01))

Question

If the velocity of projectile is not fast, it should fall down. Find out the minimal speed necessary to the moon.

In [ ]:

Stability of System ODE's

Consider 2D ODE's as follows:

\begin{eqnarray} \frac{ d x}{ d t} &=& 2 x (1-x/2 ) - x y\\ \frac{ d y}{ d t} &=& 3 y (1 -y/3) - 2 x y \end{eqnarray}
In [29]:
def f(x,y): return 2.*x*(1-x/2.) - x*y
def g(x,y): return 3.*y*(1.-y/3.) - 2.*x*y

Vector Field of Flows

It could describe how the system revolve:

In [34]:
fig1 = plt.figure(figsize=(9,6))
ax = fig1.add_subplot(1,1,1)
xmin=-0.5;  xmax=3.5;  ymin=-0.5; ymax=4.5
xmesh,ymesh=np.linspace(xmin,xmax,25),np.linspace(ymin,ymax,21)

plotwidth=5.0 # inches, horizontal distance between xmin and xmax on display
plotheight=3.9 # inches, vertical distance between ymin and ymax on display
alpha= plotwidth/(xmax-xmin);  beta = plotheight/(ymax-ymin) # scaling factors

for xp in xmesh:
  for yp in ymesh:
    vx=f(xp,yp); vy=g(xp,yp) # direction vector v=(vx,vy)
    ax.quiver(xp, yp, vx, vy, pivot='middle', headwidth=3, headlength=4,color="red")

plt.scatter((0,1,2,0),(3,1,0,0),c='b',marker='o',alpha=0.6)

ax.set_xlabel('x'); ax.set_ylabel('y')
ax.set_aspect(beta/alpha)
ax.set_xlim(xmin,xmax); ax.set_ylim(ymin,ymax)
plt.grid()                             # turn on a grid of light, dotted lines
plt.autoscale(enable=True,tight=True)  # remove extra whitespace around slopefield
plt.title("Direction Field for\n $x'= 2x(1-x/2)-xy, y'= 3y(1-y/3)-2xy$", size=16)
Out[34]:
<matplotlib.text.Text at 0x10bee4198>
/Users/cch/anaconda/lib/python3.5/site-packages/matplotlib/quiver.py:645: RuntimeWarning: divide by zero encountered in double_scalars
  length = a * (widthu_per_lenu / (self.scale * self.width))
/Users/cch/anaconda/lib/python3.5/site-packages/matplotlib/quiver.py:645: RuntimeWarning: invalid value encountered in multiply
  length = a * (widthu_per_lenu / (self.scale * self.width))

Conclusion

There are four stabilities in system;

  • the point at left top and right bottom: always stay here, quiet stable;
  • the point at left bottom: always left;
  • the point around cente: saddle point.
In [2]:
!jupyter nbconvert DE.ipynb

[NbConvertApp] Converting notebook DE.ipynb to html
[NbConvertApp] Writing 5959201 bytes to DE.html
In [ ]: